ex8

non linear regression 
with one explanatory variable  
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}
single linear regression

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)
## Initial log joint probability = -8616.18 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -41.7695    0.00198504    0.00147553      0.9751      0.9751       36    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -41.77
##    b0         7.04
##    b1         1.85
##    s          4.90
##    m0[1]     29.02
##    m0[2]     20.97
##    m0[3]     38.71
##    m0[4]     15.86
##    m0[5]     15.97
##    m0[6]     16.76
##    m0[7]     37.52
##    m0[8]     38.07
##    m0[9]     24.05
##    m0[10]    21.82
##    m0[11]    15.88
##    m0[12]    30.32
##    m0[13]    27.14
##    m0[14]    33.01
##    m0[15]    39.34
##    m0[16]     8.05
##    m0[17]    31.11
##    m0[18]    28.46
##    m0[19]    30.55
##    m0[20]    43.41
##    y0[1]     28.69
##    y0[2]     15.77
##    y0[3]     31.38
##    y0[4]     14.15
##    y0[5]     26.53
##    y0[6]     20.85
##    y0[7]     41.17
##    y0[8]     42.76
##    y0[9]     19.51
##    y0[10]    28.79
##    y0[11]    15.93
##    y0[12]    28.48
##    y0[13]    23.36
##    y0[14]    38.88
##    y0[15]    36.94
##    y0[16]    10.18
##    y0[17]    33.10
##    y0[18]    26.09
##    y0[19]    34.13
##    y0[20]    33.11
##    m1[1]      8.05
##    m1[2]     11.98
##    m1[3]     15.90
##    m1[4]     19.83
##    m1[5]     23.76
##    m1[6]     27.69
##    m1[7]     31.62
##    m1[8]     35.55
##    m1[9]     39.48
##    m1[10]    43.41
##    y1[1]      3.90
##    y1[2]     -3.07
##    y1[3]     23.56
##    y1[4]     20.25
##    y1[5]     13.94
##    y1[6]     31.82
##    y1[7]     27.56
##    y1[8]     39.82
##    y1[9]     45.52
##    y1[10]    39.72
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -41.80 -41.43 1.29 1.08 -44.34 -40.36 1.00      464      605
##    b0       7.01   6.84 2.93 2.80   2.46  12.07 1.01      297      327
##    b1       1.85   1.86 0.25 0.23   1.43   2.24 1.01      348      406
##    s        5.59   5.45 1.03 0.97   4.17   7.47 1.00     1031     1048
##    m0[1]   28.99  28.99 1.29 1.28  26.85  31.12 1.00     1605     1744
##    m0[2]   20.94  20.91 1.50 1.42  18.53  23.52 1.01      396      511
##    m0[3]   38.68  38.75 2.00 1.88  35.29  41.75 1.00      990     1113
##    m0[4]   15.83  15.78 1.95 1.85  12.77  19.14 1.01      320      363
##    m0[5]   15.94  15.89 1.94 1.83  12.90  19.22 1.01      321      367
##    m0[6]   16.73  16.68 1.86 1.74  13.81  19.92 1.01      327      356
##    m0[7]   37.50  37.59 1.88 1.79  34.32  40.40 1.00     1109     1349
##    m0[8]   38.04  38.12 1.93 1.82  34.76  41.01 1.00     1052     1235
##    m0[9]   24.02  24.02 1.33 1.26  21.87  26.27 1.00      555      609
##    m0[10]  21.79  21.78 1.44 1.38  19.40  24.23 1.01      425      557
##    m0[11]  15.85  15.80 1.95 1.85  12.79  19.15 1.01      320      363
##    m0[12]  30.29  30.31 1.34 1.32  28.07  32.46 1.00     1966     1718
##    m0[13]  27.12  27.11 1.27 1.25  24.99  29.22 1.00     1016     1329
##    m0[14]  32.99  33.04 1.49 1.45  30.52  35.31 1.00     2026     1324
##    m0[15]  39.31  39.40 2.06 1.93  35.75  42.48 1.00      935      967
##    m0[16]   8.02   7.87 2.81 2.69   3.66  12.95 1.01      297      327
##    m0[17]  31.08  31.11 1.38 1.35  28.77  33.28 1.00     2129     1643
##    m0[18]  28.43  28.42 1.28 1.26  26.31  30.54 1.00     1426     1750
##    m0[19]  30.52  30.54 1.35 1.33  28.27  32.68 1.00     2017     1718
##    m0[20]  43.38  43.45 2.51 2.34  39.10  47.25 1.00      703      744
##    y0[1]   29.02  28.98 5.79 5.60  19.65  38.30 1.00     2187     2011
##    y0[2]   20.85  20.80 5.80 5.47  11.59  30.23 1.00     1631     1518
##    y0[3]   38.75  38.67 6.08 5.84  28.99  48.64 1.00     1714     1709
##    y0[4]   15.73  15.75 5.95 5.73   6.16  25.49 1.00     1207     1483
##    y0[5]   16.05  16.11 6.03 5.85   6.26  26.11 1.00     1316     1420
##    y0[6]   16.62  16.76 5.95 5.84   6.75  26.05 1.00     1210     1243
##    y0[7]   37.60  37.52 6.06 5.85  27.80  47.60 1.00     1974     1988
##    y0[8]   38.21  38.36 5.93 5.70  28.29  47.64 1.00     1978     2007
##    y0[9]   24.14  24.12 6.06 5.61  14.42  34.05 1.00     1642     1881
##    y0[10]  21.71  21.82 5.72 5.50  12.67  30.95 1.00     1742     1671
##    y0[11]  15.77  15.82 6.01 5.76   5.99  25.68 1.00     1272     1745
##    y0[12]  30.18  30.06 5.80 5.46  20.70  39.74 1.00     2002     1597
##    y0[13]  27.09  27.05 5.90 5.53  17.38  36.61 1.00     2160     2054
##    y0[14]  32.89  32.97 5.85 5.70  23.30  42.76 1.00     1965     1900
##    y0[15]  39.41  39.37 6.13 5.79  28.97  49.30 1.00     1998     1834
##    y0[16]   8.04   7.96 6.24 6.05  -2.15  18.16 1.01      951     1250
##    y0[17]  31.15  31.10 5.88 5.65  21.69  40.58 1.00     1905     1931
##    y0[18]  28.35  28.35 5.82 5.68  18.76  37.99 1.00     1965     1894
##    y0[19]  30.63  30.55 5.68 5.43  21.39  39.89 1.00     2261     1797
##    y0[20]  43.35  43.43 6.25 5.84  33.10  53.53 1.00     1570     1729
##    m1[1]    8.02   7.87 2.81 2.69   3.66  12.95 1.01      297      327
##    m1[2]   11.95  11.87 2.36 2.28   8.16  16.05 1.01      302      336
##    m1[3]   15.88  15.82 1.94 1.84  12.82  19.17 1.01      320      363
##    m1[4]   19.81  19.78 1.59 1.49  17.30  22.55 1.01      370      437
##    m1[5]   23.73  23.73 1.34 1.28  21.55  26.00 1.00      532      620
##    m1[6]   27.66  27.65 1.27 1.26  25.54  29.75 1.00     1163     1473
##    m1[7]   31.59  31.64 1.40 1.38  29.22  33.85 1.00     2171     1588
##    m1[8]   35.52  35.61 1.70 1.62  32.70  38.17 1.00     1400     1363
##    m1[9]   39.45  39.54 2.08 1.95  35.86  42.65 1.00      924      973
##    m1[10]  43.38  43.45 2.51 2.34  39.10  47.25 1.00      703      744
##    y1[1]    8.01   7.95 6.40 6.12  -2.19  18.58 1.00     1112     1455
##    y1[2]   11.78  11.84 6.00 5.66   2.02  21.46 1.00     1420     1515
##    y1[3]   15.73  15.84 5.99 5.96   6.01  25.56 1.00     1417     1589
##    y1[4]   19.70  19.58 6.00 5.86  10.10  29.44 1.00     1540     2003
##    y1[5]   23.88  23.87 5.84 5.80  14.27  33.16 1.00     1892     1930
##    y1[6]   27.66  27.69 5.88 5.45  18.33  37.77 1.00     2068     1866
##    y1[7]   31.45  31.72 5.85 5.79  21.69  40.75 1.00     2001     1971
##    y1[8]   35.65  35.40 5.85 5.60  26.09  45.50 1.00     1879     2014
##    y1[9]   39.36  39.43 6.06 5.74  29.60  49.02 1.00     1815     1910
##    y1[10]  43.14  43.12 6.12 5.87  33.37  53.30 1.00     1480     1703

quadratic regression  
y=b0+b2(x-b1)**2  

ex8-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real b2;
  real<lower=0> s;
}
model{
  y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b2*(x[i]-b1)^2;
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b2*(x1[i]-b1)^2;
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)
## Initial log joint probability = -137.406 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       37      -9.26199     0.0001523    0.00171394           1           1       47    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -9.26
##    b0         4.16
##    b1         3.99
##    b2         0.54
##    s          0.96
##    m0[1]      5.55
##    m0[2]      5.28
##    m0[3]     15.12
##    m0[4]     14.50
##    m0[5]      9.77
##    m0[6]     11.75
##    m0[7]      5.19
##    m0[8]     10.77
##    m0[9]      5.44
##    m0[10]     7.54
##    m0[11]    17.58
##    m0[12]     8.87
##    m0[13]     6.28
##    m0[14]     4.58
##    m0[15]     6.33
##    m0[16]     8.61
##    m0[17]     6.28
##    m0[18]     9.71
##    m0[19]    10.38
##    m0[20]    14.79
##    y0[1]      7.37
##    y0[2]      6.68
##    y0[3]     15.33
##    y0[4]     13.88
##    y0[5]      9.92
##    y0[6]     12.58
##    y0[7]      4.18
##    y0[8]     11.46
##    y0[9]      3.99
##    y0[10]     9.51
##    y0[11]    18.57
##    y0[12]     8.92
##    y0[13]     5.52
##    y0[14]     5.61
##    y0[15]     6.28
##    y0[16]     8.05
##    y0[17]     7.76
##    y0[18]     9.14
##    y0[19]    11.20
##    y0[20]    13.57
##    m1[1]     11.75
##    m1[2]      8.33
##    m1[3]      5.93
##    m1[4]      4.54
##    m1[5]      4.17
##    m1[6]      4.82
##    m1[7]      6.48
##    m1[8]      9.16
##    m1[9]     12.86
##    m1[10]    17.58
##    y1[1]     12.30
##    y1[2]      8.88
##    y1[3]      3.58
##    y1[4]      3.38
##    y1[5]      3.67
##    y1[6]      4.19
##    y1[7]      4.92
##    y1[8]      8.97
##    y1[9]     13.30
##    y1[10]    17.73
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -11.46 -11.11 1.55 1.32 -14.37 -9.65 1.00      546      851
##    b0       4.15   4.16 0.42 0.41   3.42  4.82 1.01      375      405
##    b1       3.98   3.98 0.09 0.08   3.84  4.12 1.00     1580     1017
##    b2       0.54   0.54 0.04 0.04   0.48  0.61 1.01      517      620
##    s        1.14   1.10 0.22 0.20   0.85  1.54 1.01      879      889
##    m0[1]    5.55   5.55 0.39 0.38   4.91  6.17 1.01      410      455
##    m0[2]    5.27   5.28 0.37 0.37   4.65  5.87 1.01      443      533
##    m0[3]   15.11  15.11 0.50 0.48  14.29 15.92 1.00     1567     1212
##    m0[4]   14.49  14.49 0.47 0.45  13.73 15.25 1.00     1780     1270
##    m0[5]    9.76   9.75 0.47 0.45   9.02 10.52 1.00     1617     1336
##    m0[6]   11.74  11.74 0.59 0.56  10.83 12.69 1.00     1599     1337
##    m0[7]    5.19   5.19 0.40 0.39   4.52  5.83 1.01      398      464
##    m0[8]   10.76  10.77 0.34 0.33  10.22 11.31 1.00     1819     1658
##    m0[9]    5.43   5.44 0.37 0.37   4.82  6.02 1.01      458      527
##    m0[10]   7.53   7.53 0.34 0.33   6.98  8.08 1.01      550      767
##    m0[11]  17.56  17.57 0.64 0.59  16.51 18.58 1.00     1123     1086
##    m0[12]   8.86   8.85 0.42 0.40   8.19  9.58 1.00     1442     1319
##    m0[13]   6.28   6.28 0.37 0.36   5.69  6.88 1.01      442      555
##    m0[14]   4.58   4.60 0.39 0.39   3.91  5.21 1.01      394      448
##    m0[15]   6.33   6.32 0.37 0.36   5.74  6.93 1.01      445      575
##    m0[16]   8.60   8.60 0.32 0.33   8.07  9.13 1.01      744     1020
##    m0[17]   6.27   6.27 0.37 0.36   5.68  6.87 1.01      442      555
##    m0[18]   9.70   9.70 0.46 0.44   8.97 10.46 1.00     1611     1336
##    m0[19]  10.37  10.36 0.50 0.48   9.58 11.19 1.00     1661     1267
##    m0[20]  14.78  14.78 0.49 0.46  13.98 15.56 1.00     1676     1244
##    y0[1]    5.59   5.63 1.22 1.16   3.62  7.58 1.00     1337     1911
##    y0[2]    5.27   5.28 1.22 1.13   3.28  7.25 1.00     1542     1712
##    y0[3]   15.12  15.14 1.25 1.19  13.07 17.15 1.00     1980     2017
##    y0[4]   14.46  14.49 1.24 1.20  12.32 16.40 1.00     2054     1633
##    y0[5]    9.78   9.76 1.28 1.21   7.74 11.89 1.00     1927     1834
##    y0[6]   11.73  11.74 1.26 1.23   9.72 13.74 1.00     1944     1758
##    y0[7]    5.19   5.19 1.25 1.16   3.17  7.16 1.00     1548     1386
##    y0[8]   10.82  10.84 1.21 1.13   8.87 12.79 1.00     2122     1913
##    y0[9]    5.40   5.45 1.22 1.17   3.35  7.37 1.00     1556     1688
##    y0[10]   7.54   7.49 1.21 1.14   5.61  9.60 1.00     1450     1740
##    y0[11]  17.56  17.59 1.30 1.21  15.37 19.62 1.00     1928     1713
##    y0[12]   8.87   8.90 1.25 1.17   6.74 10.89 1.00     1779     1756
##    y0[13]   6.30   6.33 1.18 1.12   4.36  8.21 1.00     1982     1880
##    y0[14]   4.59   4.61 1.21 1.12   2.65  6.60 1.00     1457     1684
##    y0[15]   6.32   6.31 1.19 1.14   4.42  8.26 1.00     1238     1745
##    y0[16]   8.58   8.57 1.20 1.12   6.72 10.51 1.00     1713     1685
##    y0[17]   6.27   6.30 1.23 1.12   4.30  8.29 1.00     1707     1758
##    y0[18]   9.69   9.69 1.21 1.15   7.73 11.65 1.00     1937     1664
##    y0[19]  10.33  10.34 1.25 1.15   8.26 12.38 1.00     1825     1750
##    y0[20]  14.74  14.76 1.23 1.16  12.73 16.75 1.00     1925     1882
##    m1[1]   11.74  11.74 0.59 0.56  10.83 12.69 1.00     1599     1337
##    m1[2]    8.32   8.32 0.40 0.38   7.69  9.01 1.00     1248     1249
##    m1[3]    5.92   5.93 0.36 0.36   5.33  6.50 1.01      520      573
##    m1[4]    4.54   4.55 0.40 0.40   3.87  5.18 1.01      392      443
##    m1[5]    4.17   4.17 0.42 0.42   3.44  4.84 1.01      375      419
##    m1[6]    4.81   4.82 0.41 0.41   4.12  5.47 1.01      388      431
##    m1[7]    6.48   6.48 0.36 0.35   5.89  7.07 1.01      454      606
##    m1[8]    9.16   9.16 0.32 0.32   8.63  9.68 1.00      910     1279
##    m1[9]   12.85  12.85 0.40 0.38  12.20 13.50 1.00     2073     1664
##    m1[10]  17.56  17.57 0.64 0.59  16.51 18.58 1.00     1123     1086
##    y1[1]   11.74  11.76 1.30 1.27   9.61 13.81 1.00     2003     1870
##    y1[2]    8.28   8.27 1.21 1.12   6.30 10.28 1.00     1740     1876
##    y1[3]    5.92   5.88 1.24 1.16   4.03  8.00 1.00     1545     1858
##    y1[4]    4.56   4.54 1.22 1.16   2.57  6.62 1.00     1419     1688
##    y1[5]    4.16   4.16 1.22 1.20   2.18  6.12 1.00     1443     1649
##    y1[6]    4.83   4.81 1.23 1.21   2.76  6.87 1.00     1413     1642
##    y1[7]    6.48   6.53 1.21 1.17   4.45  8.47 1.00     1453     1719
##    y1[8]    9.15   9.15 1.22 1.19   7.12 11.18 1.00     1816     1813
##    y1[9]   12.85  12.82 1.21 1.17  10.89 14.82 1.00     1814     1855
##    y1[10]  17.56  17.56 1.33 1.31  15.39 19.67 1.00     1840     1578

both log regression  
log y=b0+b1*log x   x,y>0
y=exp b0 * x**b1

ex8-2.stan

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector<lower=0>[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~lognormal(b0+b1*log(x),s);
}
generated quantities{
  vector[N] m0;
  vector<lower=0>[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*log(x[i]);
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector<lower=0>[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*log(x1[i]);
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
             qplot(log(x),log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -5181.34 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -29.9396   0.000262506   0.000132818           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -29.94
##    b0         1.05
##    b1         2.16
##    s          1.08
##    m0[1]      5.31
##    m0[2]      5.60
##    m0[3]      1.36
##    m0[4]      5.97
##    m0[5]      2.09
##    m0[6]      5.86
##    m0[7]      2.52
##    m0[8]      5.59
##    m0[9]      6.01
##    m0[10]     3.37
##    m0[11]     4.61
##    m0[12]     4.58
##    m0[13]     5.06
##    m0[14]     5.62
##    m0[15]     5.45
##    m0[16]     5.27
##    m0[17]    -2.68
##    m0[18]     1.15
##    m0[19]     4.28
##    m0[20]     0.65
##    y0[1]     40.62
##    y0[2]    343.62
##    y0[3]      1.00
##    y0[4]     48.48
##    y0[5]      3.30
##    y0[6]     59.33
##    y0[7]     83.10
##    y0[8]   1195.92
##    y0[9]   1455.42
##    y0[10]    23.08
##    y0[11]   163.65
##    y0[12]   525.74
##    y0[13]    90.22
##    y0[14]   135.41
##    y0[15]   276.20
##    y0[16]   101.27
##    y0[17]     0.14
##    y0[18]     3.47
##    y0[19]   199.75
##    y0[20]     2.77
##    m1[1]     -2.68
##    m1[2]      1.55
##    m1[3]      2.89
##    m1[4]      3.71
##    m1[5]      4.30
##    m1[6]      4.77
##    m1[7]      5.15
##    m1[8]      5.48
##    m1[9]      5.76
##    m1[10]     6.01
##    y1[1]      0.03
##    y1[2]      7.12
##    y1[3]     21.24
##    y1[4]     60.92
##    y1[5]     78.03
##    y1[6]    412.04
##    y1[7]     70.63
##    y1[8]   2071.88
##    y1[9]    124.94
##    y1[10]   250.87
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable    mean median      sd    mad     q5     q95 rhat ess_bulk ess_tail
##    lp__    -31.50 -31.16    1.34   1.08 -34.22  -30.05 1.00      706     1024
##    b0        1.04   1.06    0.47   0.46   0.25    1.76 1.01      393      686
##    b1        2.17   2.16    0.28   0.27   1.73    2.62 1.01      478      861
##    s         1.23   1.20    0.23   0.21   0.93    1.65 1.00     1385     1425
##    m0[1]     5.31   5.31    0.32   0.31   4.78    5.82 1.01     1952     1457
##    m0[2]     5.60   5.60    0.34   0.33   5.03    6.15 1.00     1833     1394
##    m0[3]     1.35   1.37    0.43   0.42   0.61    2.03 1.01      397      713
##    m0[4]     5.98   5.98    0.37   0.37   5.36    6.57 1.00     1569     1457
##    m0[5]     2.08   2.10    0.37   0.36   1.48    2.67 1.01      412      763
##    m0[6]     5.87   5.87    0.36   0.35   5.26    6.44 1.00     1644     1374
##    m0[7]     2.51   2.52    0.33   0.33   1.98    3.05 1.01      431      855
##    m0[8]     5.60   5.60    0.34   0.33   5.03    6.14 1.00     1836     1394
##    m0[9]     6.02   6.02    0.37   0.37   5.40    6.61 1.00     1538     1435
##    m0[10]    3.37   3.37    0.28   0.27   2.91    3.83 1.01      575     1028
##    m0[11]    4.61   4.61    0.28   0.27   4.15    5.06 1.00     1448     1665
##    m0[12]    4.58   4.58    0.28   0.27   4.13    5.03 1.00     1416     1690
##    m0[13]    5.07   5.07    0.30   0.29   4.56    5.56 1.00     1888     1586
##    m0[14]    5.63   5.63    0.34   0.33   5.05    6.18 1.00     1817     1394
##    m0[15]    5.46   5.46    0.33   0.32   4.91    5.99 1.00     1917     1437
##    m0[16]    5.27   5.27    0.31   0.30   4.75    5.78 1.01     1969     1459
##    m0[17]   -2.72  -2.68    0.91   0.89  -4.24   -1.34 1.01      403      721
##    m0[18]    1.14   1.16    0.46   0.45   0.37    1.85 1.01      394      725
##    m0[19]    4.28   4.28    0.27   0.26   3.83    4.71 1.01     1067     1356
##    m0[20]    0.64   0.66    0.51   0.50  -0.24    1.41 1.01      390      671
##    y0[1]   483.09 201.70 1255.66 207.38  25.15 1601.75 1.00     1877     1932
##    y0[2]   648.86 274.92 2153.91 290.67  33.93 2209.18 1.00     2138     2044
##    y0[3]     9.08   3.86   20.12   4.02   0.44   31.41 1.00     1283     1966
##    y0[4]   974.43 382.68 2611.32 402.41  44.02 3163.10 1.00     2038     1846
##    y0[5]    18.92   8.32   65.03   8.79   0.96   61.62 1.00     1619     1973
##    y0[6]   832.19 359.83 1786.01 375.18  38.78 2957.74 1.00     1934     1712
##    y0[7]    30.65  12.70   70.65  13.96   1.40  105.57 1.00     1624     1764
##    y0[8]   806.97 260.03 9360.08 270.99  29.11 2067.98 1.00     2048     1696
##    y0[9]  1063.43 434.67 2777.23 450.57  56.50 3707.01 1.00     1868     1773
##    y0[10]   68.00  29.38  144.46  30.71   3.52  242.03 1.00     1677     1970
##    y0[11]  238.16  99.98  731.96 103.69  12.46  757.54 1.00     1978     1856
##    y0[12]  229.47  94.96  482.03  95.30  12.06  821.65 1.00     2101     1799
##    y0[13]  366.09 160.51  806.35 168.53  18.40 1382.51 1.00     2035     1945
##    y0[14]  702.27 269.73 2291.80 278.49  32.49 2346.89 1.00     1837     1647
##    y0[15]  531.43 222.82 1377.34 233.42  25.85 1760.23 1.00     2065     1891
##    y0[16]  499.38 192.34 1380.56 198.56  24.29 1770.73 1.00     2169     1927
##    y0[17]    0.24   0.07    0.97   0.08   0.01    0.85 1.01      814     1322
##    y0[18]    7.58   3.25   17.59   3.53   0.35   27.53 1.00     1473     1831
##    y0[19]  160.03  70.96  283.27  74.62   7.78  571.37 1.00     1623     1875
##    y0[20]    4.50   1.84   10.33   2.05   0.19   15.70 1.00     1720     1593
##    m1[1]    -2.72  -2.68    0.91   0.89  -4.24   -1.34 1.01      403      721
##    m1[2]     1.54   1.55    0.42   0.41   0.84    2.19 1.01      401      708
##    m1[3]     2.88   2.89    0.30   0.30   2.39    3.39 1.01      467      992
##    m1[4]     3.71   3.71    0.27   0.26   3.26    4.15 1.01      702     1101
##    m1[5]     4.30   4.30    0.27   0.26   3.86    4.74 1.01     1088     1425
##    m1[6]     4.77   4.77    0.29   0.28   4.30    5.24 1.00     1621     1619
##    m1[7]     5.16   5.16    0.31   0.29   4.64    5.66 1.00     1930     1587
##    m1[8]     5.48   5.48    0.33   0.32   4.93    6.02 1.00     1905     1437
##    m1[9]     5.77   5.77    0.35   0.35   5.17    6.33 1.00     1718     1375
##    m1[10]    6.02   6.02    0.37   0.37   5.40    6.61 1.00     1538     1435
##    y1[1]     0.23   0.06    0.75   0.08   0.01    0.75 1.01      977     1254
##    y1[2]    14.27   4.60  136.39   4.64   0.55   40.32 1.00     1531     1935
##    y1[3]    46.27  17.50  125.87  18.73   2.01  162.26 1.00     1652     1892
##    y1[4]    97.46  39.66  238.91  41.78   4.62  305.60 1.00     2085     1918
##    y1[5]   176.58  77.52  418.60  81.21   8.36  596.34 1.00     1951     1972
##    y1[6]   307.66 110.52  982.07 117.52  13.80  968.40 1.00     2087     1988
##    y1[7]   412.31 157.17 1033.51 165.26  20.78 1434.79 1.00     1995     1845
##    y1[8]   600.40 253.31 1467.35 268.70  32.10 2043.06 1.00     1944     1890
##    y1[9]   805.42 322.38 1953.98 325.43  41.24 2824.73 1.00     1996     1914
##    y1[10]  986.18 394.15 2306.50 414.51  50.62 3518.26 1.00     2016     1848

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

qplot(log(x),log(y),col=I('red'))+
  geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=m),xy,col='black')

ex8-3

exponential increasing/decreasing  

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  
log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)  
b1>0 x->Infinity,y->Infinity  
b1<0 x->Infinity,y->+0  
n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  

ex8-3-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0*exp(b1*x),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*exp(b1*x[i]);
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*exp(b1*x1[i]);
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)
## Initial log joint probability = -24.6144 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19       10.8609   4.19475e-05    0.00133438           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      10.86
##    b0        15.08
##    b1        -2.34
##    s          0.35
##    m0[1]      0.00
##    m0[2]      0.02
##    m0[3]      3.50
##    m0[4]      0.00
##    m0[5]      0.00
##    m0[6]      0.03
##    m0[7]      0.00
##    m0[8]      0.37
##    m0[9]      0.00
##    m0[10]     1.08
##    m0[11]     0.11
##    m0[12]     0.00
##    m0[13]     0.12
##    m0[14]     1.54
##    m0[15]     0.69
##    m0[16]     1.89
##    m0[17]     1.92
##    m0[18]     0.00
##    m0[19]     0.68
##    m0[20]     0.09
##    y0[1]     -0.44
##    y0[2]      0.67
##    y0[3]      3.91
##    y0[4]      0.09
##    y0[5]      0.01
##    y0[6]      0.09
##    y0[7]      0.07
##    y0[8]      0.49
##    y0[9]     -0.04
##    y0[10]     0.58
##    y0[11]     0.22
##    y0[12]    -0.15
##    y0[13]     0.18
##    y0[14]     1.26
##    y0[15]     1.58
##    y0[16]     1.95
##    y0[17]     1.97
##    y0[18]     0.13
##    y0[19]     0.96
##    y0[20]    -0.03
##    m1[1]      3.50
##    m1[2]      1.28
##    m1[3]      0.47
##    m1[4]      0.17
##    m1[5]      0.06
##    m1[6]      0.02
##    m1[7]      0.01
##    m1[8]      0.00
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]      3.74
##    y1[2]      1.40
##    y1[3]      0.11
##    y1[4]      0.13
##    y1[5]      0.17
##    y1[6]     -0.01
##    y1[7]      0.07
##    y1[8]      0.16
##    y1[9]     -0.64
##    y1[10]    -0.03
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   12.42  12.77 1.40 1.13  9.79 13.93 1.00      477      436
##    b0     19.45  17.18 9.28 5.97 10.39 35.30 1.02      287      234
##    b1     -2.58  -2.49 0.50 0.44 -3.47 -1.90 1.01      285      255
##    s       0.41   0.39 0.08 0.07  0.30  0.55 1.00      604      731
##    m0[1]   0.00   0.00 0.00 0.00  0.00  0.01 1.01      289      273
##    m0[2]   0.02   0.02 0.02 0.01  0.00  0.06 1.01      291      302
##    m0[3]   3.60   3.58 0.40 0.39  2.97  4.27 1.00      489      459
##    m0[4]   0.00   0.00 0.00 0.00  0.00  0.01 1.01      289      285
##    m0[5]   0.00   0.00 0.00 0.00  0.00  0.01 1.01      289      285
##    m0[6]   0.03   0.02 0.02 0.02  0.00  0.07 1.01      291      302
##    m0[7]   0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    m0[8]   0.32   0.32 0.12 0.12  0.14  0.53 1.01      305      325
##    m0[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    m0[10]  1.00   1.01 0.19 0.19  0.68  1.30 1.01      360      335
##    m0[11]  0.10   0.09 0.06 0.05  0.02  0.20 1.01      295      294
##    m0[12]  0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    m0[13]  0.10   0.09 0.06 0.05  0.03  0.21 1.01      295      294
##    m0[14]  1.46   1.47 0.20 0.19  1.13  1.78 1.01      472      388
##    m0[15]  0.62   0.62 0.17 0.16  0.35  0.89 1.01      321      339
##    m0[16]  1.83   1.83 0.20 0.18  1.49  2.14 1.00      694      598
##    m0[17]  1.85   1.86 0.20 0.18  1.52  2.17 1.00      722      617
##    m0[18]  0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    m0[19]  0.61   0.61 0.16 0.16  0.34  0.88 1.01      321      339
##    m0[20]  0.08   0.07 0.05 0.04  0.02  0.17 1.01      294      292
##    y0[1]   0.00   0.00 0.40 0.39 -0.66  0.65 1.00     1735     1797
##    y0[2]   0.00  -0.01 0.43 0.41 -0.69  0.70 1.00     1870     1893
##    y0[3]   3.61   3.60 0.58 0.55  2.70  4.54 1.00      895      972
##    y0[4]  -0.01   0.00 0.42 0.40 -0.68  0.66 1.00     2075     1995
##    y0[5]  -0.01  -0.01 0.41 0.39 -0.68  0.66 1.00     1998     2012
##    y0[6]   0.04   0.04 0.41 0.38 -0.62  0.72 1.00     1920     2053
##    y0[7]  -0.01  -0.01 0.41 0.39 -0.69  0.66 1.00     2064     2052
##    y0[8]   0.33   0.34 0.41 0.40 -0.36  0.96 1.00     1646     1809
##    y0[9]  -0.01   0.00 0.41 0.38 -0.67  0.66 1.00     1992     1591
##    y0[10]  0.97   0.98 0.45 0.42  0.22  1.72 1.00     1148     1248
##    y0[11]  0.10   0.11 0.41 0.39 -0.59  0.79 1.00     1886     1887
##    y0[12]  0.00   0.00 0.42 0.39 -0.68  0.68 1.00     2052     1940
##    y0[13]  0.12   0.12 0.42 0.40 -0.56  0.83 1.00     1918     1766
##    y0[14]  1.46   1.46 0.47 0.43  0.67  2.22 1.00     1664     1895
##    y0[15]  0.61   0.61 0.43 0.42 -0.09  1.32 1.00     1410     1696
##    y0[16]  1.82   1.84 0.47 0.44  1.03  2.57 1.00     1433     1785
##    y0[17]  1.85   1.86 0.45 0.45  1.08  2.57 1.00     1784     1512
##    y0[18]  0.01   0.01 0.42 0.40 -0.67  0.72 1.00     2026     1893
##    y0[19]  0.62   0.62 0.45 0.43 -0.13  1.33 1.00     1456     1280
##    y0[20]  0.08   0.08 0.41 0.40 -0.59  0.75 1.00     2209     1837
##    m1[1]   3.60   3.58 0.40 0.39  2.97  4.27 1.00      489      459
##    m1[2]   1.20   1.20 0.19 0.19  0.87  1.50 1.01      394      348
##    m1[3]   0.42   0.41 0.14 0.14  0.20  0.65 1.01      309      340
##    m1[4]   0.15   0.14 0.08 0.07  0.05  0.29 1.01      297      303
##    m1[5]   0.06   0.05 0.04 0.03  0.01  0.13 1.01      293      292
##    m1[6]   0.02   0.02 0.02 0.01  0.00  0.06 1.01      291      302
##    m1[7]   0.01   0.01 0.01 0.01  0.00  0.02 1.01      290      285
##    m1[8]   0.00   0.00 0.01 0.00  0.00  0.01 1.01      289      285
##    m1[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    m1[10]  0.00   0.00 0.00 0.00  0.00  0.00 1.01      288      274
##    y1[1]   3.59   3.60 0.57 0.57  2.67  4.53 1.00      864      908
##    y1[2]   1.20   1.21 0.47 0.45  0.43  1.94 1.00     1374     1336
##    y1[3]   0.41   0.42 0.44 0.41 -0.33  1.13 1.00     1519     1678
##    y1[4]   0.15   0.15 0.42 0.39 -0.52  0.83 1.00     1974     2055
##    y1[5]   0.06   0.05 0.42 0.41 -0.64  0.74 1.00     1975     1562
##    y1[6]   0.03   0.04 0.42 0.42 -0.68  0.69 1.00     1945     1811
##    y1[7]   0.01   0.00 0.42 0.40 -0.69  0.70 1.00     2025     1625
##    y1[8]   0.03   0.04 0.42 0.40 -0.65  0.72 1.00     1938     1840
##    y1[9]   0.02   0.02 0.41 0.40 -0.61  0.72 1.00     2338     1857
##    y1[10]  0.00   0.02 0.41 0.40 -0.68  0.64 1.00     1934     1895

log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)

ex8-3-2.stan

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real<lower=-10,upper=10> b1;
  real<lower=0> s;
}
model{
  y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=log(b0)+b1*x[i];
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=log(b0)+b1*x1[i];
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -14383.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12      -79.5531    0.00216098     0.0003441      0.9588      0.9588       31    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable                estimate
##    lp__                    -79.55
##    b0              17121900000.00
##    b1                      -10.00
##    s                        12.92
##    m0[1]                    16.93
##    m0[2]                   -25.28
##    m0[3]                    10.32
##    m0[4]                   -17.48
##    m0[5]                   -25.54
##    m0[6]                     8.44
##    m0[7]                    23.15
##    m0[8]                   -18.20
##    m0[9]                   -21.39
##    m0[10]                  -16.04
##    m0[11]                  -11.34
##    m0[12]                    9.98
##    m0[13]                  -18.46
##    m0[14]                    6.64
##    m0[15]                  -20.83
##    m0[16]                    1.48
##    m0[17]                  -10.69
##    m0[18]                   18.36
##    m0[19]                   18.71
##    m0[20]                    4.04
##    y0[1]                     1.37
##    y0[2]                     2.88
##    y0[3]          168823000000.00
##    y0[4]                     0.00
##    y0[5]                     0.00
##    y0[6]               6574500.00
##    y0[7]                     0.00
##    y0[8]                     0.00
##    y0[9]                     0.00
##    y0[10]                    0.00
##    y0[11]                    0.02
##    y0[12]             30023400.00
##    y0[13]                    0.09
##    y0[14]               286488.00
##    y0[15]                    0.00
##    y0[16]                    8.80
##    y0[17]                    0.00
##    y0[18]          14414600000.00
##    y0[19]          13937800000.00
##    y0[20]                    0.00
##    m1[1]                    23.15
##    m1[2]                    17.74
##    m1[3]                    12.33
##    m1[4]                     6.92
##    m1[5]                     1.51
##    m1[6]                    -3.90
##    m1[7]                    -9.31
##    m1[8]                   -14.72
##    m1[9]                   -20.13
##    m1[10]                  -25.54
##    y1[1]              67628300.00
##    y1[2]                 13202.10
##    y1[3]             792279000.00
##    y1[4]  10221900000000000000.00
##    y1[5]                     0.00
##    y1[6]                     0.03
##    y1[7]                     0.00
##    y1[8]                     2.50
##    y1[9]                     0.00
##    y1[10]                    0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   -7.40  -7.08 1.26 1.04 -9.78 -6.05 1.01      654      829
##    b0      9.36   9.09 1.72 1.55  6.96 12.49 1.01      650      473
##    b1     -2.07  -2.07 0.06 0.05 -2.17 -1.98 1.00      831      700
##    s       0.41   0.40 0.07 0.06  0.31  0.55 1.01      686      664
##    m0[1]   0.85   0.84 0.15 0.14  0.62  1.10 1.01      655      463
##    m0[2]  -7.89  -7.88 0.16 0.15 -8.14 -7.63 1.00     1811     1327
##    m0[3]  -0.52  -0.53 0.12 0.12 -0.71 -0.31 1.00      699      525
##    m0[4]  -6.27  -6.27 0.13 0.12 -6.48 -6.07 1.00     2134     1590
##    m0[5]  -7.94  -7.94 0.16 0.16 -8.20 -7.68 1.00     1800     1326
##    m0[6]  -0.91  -0.91 0.11 0.11 -1.09 -0.71 1.00      726      597
##    m0[7]   2.13   2.12 0.17 0.17  1.86  2.43 1.01      650      473
##    m0[8]  -6.42  -6.42 0.13 0.13 -6.63 -6.21 1.00     2111     1445
##    m0[9]  -7.08  -7.08 0.14 0.14 -7.31 -6.86 1.00     1980     1425
##    m0[10] -5.97  -5.97 0.12 0.12 -6.17 -5.78 1.00     2165     1374
##    m0[11] -5.00  -5.00 0.11 0.11 -5.17 -4.83 1.00     2153     1508
##    m0[12] -0.59  -0.60 0.12 0.12 -0.78 -0.38 1.00      704      525
##    m0[13] -6.47  -6.47 0.13 0.13 -6.69 -6.27 1.00     2102     1444
##    m0[14] -1.28  -1.29 0.11 0.10 -1.45 -1.09 1.00      758      586
##    m0[15] -6.96  -6.96 0.14 0.14 -7.19 -6.75 1.00     2004     1403
##    m0[16] -2.35  -2.35 0.10 0.09 -2.51 -2.18 1.00      931      802
##    m0[17] -4.87  -4.87 0.10 0.10 -5.04 -4.70 1.00     2117     1555
##    m0[18]  1.14   1.14 0.15 0.15  0.90  1.41 1.01      652      475
##    m0[19]  1.22   1.21 0.15 0.15  0.97  1.48 1.01      651      475
##    m0[20] -1.82  -1.82 0.10 0.10 -1.98 -1.65 1.00      824      676
##    y0[1]   2.62   2.34 1.24 0.96  1.19  5.03 1.00     1710     1618
##    y0[2]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     2073     1550
##    y0[3]   0.65   0.59 0.31 0.25  0.28  1.21 1.00     1690     1846
##    y0[4]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     2103     1933
##    y0[5]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     2050     1884
##    y0[6]   0.45   0.40 0.21 0.16  0.20  0.85 1.00     1988     1970
##    y0[7]   9.43   8.41 4.87 3.62  4.04 18.50 1.00     1352     1470
##    y0[8]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     1865     1955
##    y0[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     2062     1945
##    y0[10]  0.00   0.00 0.00 0.00  0.00  0.01 1.00     1854     1714
##    y0[11]  0.01   0.01 0.00 0.00  0.00  0.01 1.00     1960     1721
##    y0[12]  0.60   0.55 0.27 0.22  0.28  1.07 1.00     1858     1855
##    y0[13]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1962     1841
##    y0[14]  0.30   0.27 0.13 0.11  0.14  0.54 1.00     1806     1893
##    y0[15]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1992     1969
##    y0[16]  0.11   0.10 0.05 0.04  0.05  0.19 1.00     1751     1839
##    y0[17]  0.01   0.01 0.00 0.00  0.00  0.02 1.00     2095     2011
##    y0[18]  3.49   3.13 1.79 1.32  1.42  6.65 1.00     1480     1608
##    y0[19]  3.70   3.36 1.80 1.41  1.57  7.10 1.00     1420     1922
##    y0[20]  0.18   0.16 0.08 0.06  0.08  0.33 1.00     1764     1567
##    m1[1]   2.13   2.12 0.17 0.17  1.86  2.43 1.01      650      473
##    m1[2]   1.02   1.01 0.15 0.15  0.78  1.27 1.01      653      463
##    m1[3]  -0.10  -0.11 0.13 0.12 -0.31  0.11 1.00      679      500
##    m1[4]  -1.22  -1.23 0.11 0.11 -1.39 -1.04 1.00      753      580
##    m1[5]  -2.34  -2.34 0.10 0.09 -2.50 -2.18 1.00      929      802
##    m1[6]  -3.46  -3.46 0.09 0.09 -3.62 -3.31 1.00     1398     1537
##    m1[7]  -4.58  -4.58 0.10 0.10 -4.74 -4.42 1.00     2020     1718
##    m1[8]  -5.70  -5.70 0.12 0.11 -5.89 -5.52 1.00     2187     1439
##    m1[9]  -6.82  -6.82 0.14 0.13 -7.04 -6.61 1.00     2035     1362
##    m1[10] -7.94  -7.94 0.16 0.16 -8.20 -7.68 1.00     1800     1326
##    y1[1]   9.21   8.34 4.59 3.53  3.88 16.99 1.00     1495     1939
##    y1[2]   3.03   2.73 1.53 1.17  1.35  5.67 1.00     1864     1709
##    y1[3]   0.99   0.90 0.47 0.38  0.43  1.82 1.00     1795     1820
##    y1[4]   0.32   0.29 0.15 0.12  0.15  0.60 1.01     1791     1910
##    y1[5]   0.11   0.10 0.05 0.04  0.05  0.20 1.00     1854     1970
##    y1[6]   0.03   0.03 0.02 0.01  0.02  0.06 1.00     1597     1949
##    y1[7]   0.01   0.01 0.01 0.00  0.01  0.02 1.00     1821     1651
##    y1[8]   0.00   0.00 0.00 0.00  0.00  0.01 1.00     1911     1919
##    y1[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.00     1757     2012
##    y1[10]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1763     1889

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,log(y),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex8-4

growth curve  

y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)  
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)
data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=10> b1;
  real<lower=0> s;
}
model{
  y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(1-exp(-b1*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(1-exp(-b1*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)
## Initial log joint probability = -222558 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -51.0573   0.000272808   0.000222331      0.9729      0.9729       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -51.06
##    b0         0.00
##    b1         0.00
##    s          7.79
##    m0[1]      0.00
##    m0[2]      0.00
##    m0[3]      0.00
##    m0[4]      0.00
##    m0[5]      0.00
##    m0[6]      0.00
##    m0[7]      0.00
##    m0[8]      0.00
##    m0[9]      0.00
##    m0[10]     0.00
##    m0[11]     0.00
##    m0[12]     0.00
##    m0[13]     0.00
##    m0[14]     0.00
##    m0[15]     0.00
##    m0[16]     0.00
##    m0[17]     0.00
##    m0[18]     0.00
##    m0[19]     0.00
##    m0[20]     0.00
##    y0[1]    -13.21
##    y0[2]     -7.29
##    y0[3]     -8.49
##    y0[4]     -2.52
##    y0[5]     14.14
##    y0[6]     14.25
##    y0[7]    -10.18
##    y0[8]     -6.30
##    y0[9]     -1.78
##    y0[10]     5.56
##    y0[11]    -9.93
##    y0[12]   -13.70
##    y0[13]     3.12
##    y0[14]     1.31
##    y0[15]     4.51
##    y0[16]     7.46
##    y0[17]    -5.70
##    y0[18]    -4.72
##    y0[19]     3.04
##    y0[20]     5.99
##    m1[1]      0.00
##    m1[2]      0.00
##    m1[3]      0.00
##    m1[4]      0.00
##    m1[5]      0.00
##    m1[6]      0.00
##    m1[7]      0.00
##    m1[8]      0.00
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]    -10.61
##    y1[2]      2.48
##    y1[3]      0.98
##    y1[4]     -1.02
##    y1[5]     -6.68
##    y1[6]    -12.56
##    y1[7]     11.68
##    y1[8]      4.19
##    y1[9]      3.50
##    y1[10]     3.16
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   -5.43  -5.09 1.31 1.07 -8.10 -4.01 1.01      688      795
##    b0     10.53  10.52 0.47 0.43  9.82 11.31 1.00      923     1068
##    b1      0.44   0.44 0.06 0.05  0.36  0.54 1.00      762      768
##    s       0.89   0.87 0.16 0.14  0.67  1.17 1.01      570      499
##    m0[1]   4.51   4.50 0.29 0.27  4.05  5.01 1.00      886      891
##    m0[2]  10.27  10.27 0.36 0.36  9.69 10.85 1.00     1107     1174
##    m0[3]   5.96   5.95 0.31 0.30  5.45  6.49 1.00      953      919
##    m0[4]   1.84   1.83 0.16 0.14  1.60  2.10 1.00      841      863
##    m0[5]   9.87   9.87 0.28 0.28  9.40 10.31 1.00     1689     1428
##    m0[6]   8.17   8.17 0.26 0.26  7.73  8.60 1.00     1468     1153
##    m0[7]   8.39   8.39 0.26 0.25  7.97  8.79 1.00     1628     1248
##    m0[8]   4.72   4.71 0.30 0.28  4.25  5.23 1.00      892      891
##    m0[9]   9.58   9.59 0.25 0.24  9.18  9.99 1.00     2130     1649
##    m0[10] 10.27  10.27 0.36 0.36  9.68 10.85 1.00     1107     1174
##    m0[11]  8.02   8.02 0.27 0.27  7.57  8.45 1.00     1375     1126
##    m0[12]  1.29   1.28 0.11 0.10  1.12  1.48 1.00      836      840
##    m0[13]  1.19   1.19 0.11 0.10  1.03  1.37 1.00      835      840
##    m0[14] 10.23  10.23 0.35 0.35  9.66 10.79 1.00     1149     1189
##    m0[15]  9.48   9.48 0.25 0.24  9.08  9.88 1.00     2259     1568
##    m0[16]  9.96   9.97 0.29 0.29  9.48 10.43 1.00     1490     1451
##    m0[17]  9.74   9.75 0.27 0.26  9.31 10.17 1.00     1886     1469
##    m0[18]  8.82   8.82 0.24 0.23  8.42  9.21 1.00     2051     1648
##    m0[19]  5.24   5.23 0.31 0.29  4.75  5.77 1.00      911      905
##    m0[20]  4.42   4.42 0.29 0.27  3.97  4.92 1.00      883      891
##    y0[1]   4.52   4.49 0.93 0.88  3.07  5.98 1.00     1661     1865
##    y0[2]  10.28  10.27 0.96 0.94  8.72 11.86 1.00     1945     1711
##    y0[3]   5.95   5.96 0.96 0.91  4.43  7.52 1.00     1807     1768
##    y0[4]   1.85   1.83 0.92 0.89  0.38  3.36 1.00     1827     1802
##    y0[5]   9.86   9.86 0.96 0.96  8.29 11.42 1.00     1888     1789
##    y0[6]   8.15   8.15 0.92 0.85  6.62  9.65 1.00     1903     1931
##    y0[7]   8.42   8.41 0.95 0.90  6.88  9.97 1.00     1956     1717
##    y0[8]   4.71   4.72 0.95 0.89  3.15  6.27 1.00     1886     1648
##    y0[9]   9.59   9.58 0.95 0.91  8.03 11.07 1.00     1986     1953
##    y0[10] 10.27  10.26 0.98 0.94  8.68 11.88 1.00     1809     1711
##    y0[11]  7.99   8.01 0.95 0.91  6.39  9.57 1.00     1859     1621
##    y0[12]  1.30   1.29 0.92 0.86 -0.21  2.78 1.00     1827     1742
##    y0[13]  1.21   1.25 0.91 0.91 -0.36  2.68 1.00     1967     1721
##    y0[14] 10.23  10.21 0.99 0.96  8.64 11.86 1.00     1937     1736
##    y0[15]  9.48   9.47 0.94 0.91  7.94 11.05 1.00     1949     1837
##    y0[16]  9.99   9.95 0.93 0.90  8.47 11.51 1.00     1782     1729
##    y0[17]  9.73   9.72 0.98 0.94  8.16 11.36 1.00     1926     1809
##    y0[18]  8.81   8.79 0.95 0.91  7.25 10.35 1.00     1949     1708
##    y0[19]  5.24   5.26 0.95 0.87  3.63  6.81 1.00     1780     1774
##    y0[20]  4.38   4.39 0.97 0.94  2.83  5.91 1.00     1879     1486
##    m1[1]   1.19   1.19 0.11 0.10  1.03  1.37 1.00      835      840
##    m1[2]   4.33   4.32 0.29 0.27  3.88  4.82 1.00      881      874
##    m1[3]   6.40   6.40 0.31 0.30  5.89  6.92 1.00      989     1039
##    m1[4]   7.77   7.78 0.28 0.27  7.31  8.23 1.00     1262     1005
##    m1[5]   8.68   8.69 0.25 0.24  8.27  9.08 1.00     1908     1583
##    m1[6]   9.29   9.30 0.24 0.24  8.89  9.69 1.00     2367     1495
##    m1[7]   9.70   9.70 0.26 0.26  9.27 10.12 1.00     1959     1465
##    m1[8]   9.97   9.98 0.30 0.29  9.49 10.44 1.00     1483     1429
##    m1[9]  10.15  10.15 0.33 0.33  9.61 10.69 1.00     1230     1259
##    m1[10] 10.27  10.27 0.36 0.36  9.69 10.85 1.00     1107     1174
##    y1[1]   1.19   1.17 0.91 0.87 -0.30  2.70 1.00     1758     2090
##    y1[2]   4.32   4.33 0.95 0.93  2.78  5.85 1.00     1813     1812
##    y1[3]   6.41   6.41 0.95 0.89  4.84  7.89 1.00     1926     2013
##    y1[4]   7.76   7.76 0.94 0.91  6.24  9.33 1.00     1934     1947
##    y1[5]   8.69   8.68 0.95 0.91  7.13 10.23 1.00     1899     1957
##    y1[6]   9.27   9.25 0.96 0.91  7.76 10.91 1.00     1991     1825
##    y1[7]   9.68   9.66 0.93 0.91  8.20 11.19 1.00     1945     1802
##    y1[8]   9.98  10.01 0.98 0.91  8.38 11.60 1.00     2004     1916
##    y1[9]  10.13  10.13 0.97 0.95  8.56 11.68 1.00     1935     1918
##    y1[10] 10.28  10.27 1.04 1.03  8.61 11.99 1.00     1730     1932

ex8-5

sigmoid curve

y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)
data{
  int N;
  vector[N] x;
  int Ym;
  array[N] int y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=100> b1;
}
model{
  y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
  array[N] int y0;
  for(i in 1:N){
    y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
  }
  array[N1] int y1;
  for(i in 1:N1){
    y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
  }
}
n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-5.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "y0"   "y1"
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -37.93 -37.61 0.99 0.69 -39.92 -37.00 1.00      591      734
##    b0       3.99   3.99 0.16 0.16   3.73   4.25 1.00     2032     1721
##    b1       1.86   1.83 0.31 0.29   1.40   2.40 1.01      390      357
##    y0[1]    0.05   0.00 0.22 0.00   0.00   0.00 1.00     1523     1503
##    y0[2]    9.97  10.00 0.17 0.00  10.00  10.00 1.00     1629       NA
##    y0[3]    8.17   8.00 1.32 1.48   6.00  10.00 1.00     1964       NA
##    y0[4]    6.64   7.00 1.59 1.48   4.00   9.00 1.00     1981     1811
##    y0[5]    9.24   9.00 0.87 1.48   8.00  10.00 1.00     1396       NA
##    y0[6]    9.98  10.00 0.16 0.00  10.00  10.00 1.00     1998       NA
##    y0[7]    0.28   0.00 0.54 0.00   0.00   1.00 1.00     1792     1756
##    y0[8]    5.63   6.00 1.69 1.48   3.00   8.00 1.00     1862     1766
##    y0[9]    0.13   0.00 0.37 0.00   0.00   1.00 1.00     1839     1898
##    y0[10]   0.10   0.00 0.33 0.00   0.00   1.00 1.00     2070     2037
##    y0[11]   9.97  10.00 0.17 0.00  10.00  10.00 1.00     2122       NA
##    y0[12]   9.73  10.00 0.53 0.00   9.00  10.00 1.00     1976       NA
##    y0[13]   0.47   0.00 0.72 0.00   0.00   2.00 1.00     1256     1383
##    y0[14]   4.41   4.00 1.69 1.48   2.00   7.00 1.00     2101     2033
##    y0[15]   9.97  10.00 0.16 0.00  10.00  10.00 1.00     1686       NA
##    y0[16]   0.20   0.00 0.47 0.00   0.00   1.00 1.00     1760     1788
##    y0[17]   9.99  10.00 0.12 0.00  10.00  10.00 1.00     2073       NA
##    y0[18]   9.95  10.00 0.23 0.00   9.00  10.00 1.00     1896       NA
##    y0[19]   8.97   9.00 1.04 1.48   7.00  10.00 1.00     1715       NA
##    y0[20]   0.12   0.00 0.36 0.00   0.00   1.00 1.00     1720     1690
##    y1[1]    0.05   0.00 0.23 0.00   0.00   0.00 1.00     1552     1560
##    y1[2]    0.17   0.00 0.44 0.00   0.00   1.00 1.00     1546     1609
##    y1[3]    0.61   0.00 0.83 0.00   0.00   2.00 1.00     1451     1528
##    y1[4]    1.99   2.00 1.38 1.48   0.00   5.00 1.00     1650     1899
##    y1[5]    5.07   5.00 1.72 1.48   2.00   8.00 1.00     1674     1753
##    y1[6]    8.08   8.00 1.36 1.48   6.00  10.00 1.00     1893       NA
##    y1[7]    9.42  10.00 0.77 0.00   8.00  10.00 1.00     1534       NA
##    y1[8]    9.84  10.00 0.41 0.00   9.00  10.00 1.00     1696       NA
##    y1[9]    9.95  10.00 0.22 0.00  10.00  10.00 1.00     1862       NA
##    y1[10]   9.99  10.00 0.12 0.00  10.00  10.00 1.00     1795       NA

y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=ymed),xy,col='black')

ex8-6

up and down

y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))
data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=200> b0;
  real<lower=0,upper=1> b1;
  real<lower=0,upper=1> b2;
  real<lower=0,upper=10> s;
}
model{
  y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-6.stan')

fn(mdl,data)
## Initial log joint probability = -454.723 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       86      -6.55155    4.3167e-06    0.00416671      0.7339      0.7339      130    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -6.55
##    b0       104.18
##    b1         0.03
##    b2         0.19
##    s          0.84
##    m0[1]     52.85
##    m0[2]     30.98
##    m0[3]     57.29
##    m0[4]     40.25
##    m0[5]     57.13
##    m0[6]     36.59
##    m0[7]     26.90
##    m0[8]     51.95
##    m0[9]     44.02
##    m0[10]    25.50
##    m0[11]    36.81
##    m0[12]    52.08
##    m0[13]    39.74
##    m0[14]    35.14
##    m0[15]    34.68
##    m0[16]    55.62
##    m0[17]    57.06
##    m0[18]    58.09
##    m0[19]    59.44
##    m0[20]    44.26
##    y0[1]     54.53
##    y0[2]     31.69
##    y0[3]     56.37
##    y0[4]     40.12
##    y0[5]     56.28
##    y0[6]     37.88
##    y0[7]     27.08
##    y0[8]     50.76
##    y0[9]     44.39
##    y0[10]    24.05
##    y0[11]    35.85
##    y0[12]    53.21
##    y0[13]    39.56
##    y0[14]    34.84
##    y0[15]    34.50
##    y0[16]    55.36
##    y0[17]    57.63
##    y0[18]    57.86
##    y0[19]    59.12
##    y0[20]    44.67
##    m1[1]     52.08
##    m1[2]     60.07
##    m1[3]     58.99
##    m1[4]     54.45
##    m1[5]     48.86
##    m1[6]     43.26
##    m1[7]     38.04
##    m1[8]     33.34
##    m1[9]     29.17
##    m1[10]    25.50
##    y1[1]     51.66
##    y1[2]     60.22
##    y1[3]     59.56
##    y1[4]     53.04
##    y1[5]     48.30
##    y1[6]     42.19
##    y1[7]     38.28
##    y1[8]     32.44
##    y1[9]     28.92
##    y1[10]    23.97
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -10.87 -10.51 1.77 1.57 -14.43  -8.74 1.04       84      311
##    b0     105.12 105.04 4.60 4.65  97.71 113.16 1.10       27      120
##    b1       0.03   0.03 0.00 0.00   0.03   0.03 1.11       26      151
##    b2       0.18   0.18 0.01 0.01   0.17   0.20 1.09       29       52
##    s        1.04   1.02 0.20 0.17   0.76   1.43 1.14       23       82
##    m0[1]   52.89  52.88 0.44 0.43  52.18  53.61 1.05       65      658
##    m0[2]   30.93  30.92 0.40 0.39  30.27  31.60 1.06       52      497
##    m0[3]   57.28  57.26 0.51 0.47  56.45  58.14 1.01      517      342
##    m0[4]   40.24  40.24 0.29 0.27  39.76  40.72 1.00     1200     1034
##    m0[5]   57.12  57.10 0.51 0.48  56.28  57.99 1.01      495      318
##    m0[6]   36.56  36.56 0.31 0.29  36.03  37.06 1.02      944      713
##    m0[7]   26.83  26.82 0.48 0.47  26.08  27.64 1.08       36      617
##    m0[8]   51.99  51.98 0.43 0.42  51.31  52.68 1.04       69      570
##    m0[9]   44.03  44.03 0.32 0.30  43.53  44.56 1.01      658      744
##    m0[10]  25.42  25.42 0.50 0.49  24.63  26.29 1.08       34      658
##    m0[11]  36.79  36.79 0.31 0.29  36.26  37.28 1.02      999      724
##    m0[12]  52.06  52.04 0.63 0.58  51.03  53.17 1.03      120      196
##    m0[13]  39.73  39.73 0.29 0.27  39.25  40.22 1.00     1232      901
##    m0[14]  35.11  35.10 0.33 0.31  34.56  35.63 1.03      188      674
##    m0[15]  34.65  34.65 0.34 0.31  34.10  35.18 1.03      144      674
##    m0[16]  55.66  55.66 0.46 0.46  54.93  56.42 1.05       62      589
##    m0[17]  57.06  57.04 0.51 0.48  56.21  57.93 1.01      488      318
##    m0[18]  58.13  58.13 0.46 0.45  57.40  58.86 1.05       71      914
##    m0[19]  59.48  59.48 0.44 0.43  58.77  60.20 1.04       96     1003
##    m0[20]  44.27  44.27 0.32 0.31  43.76  44.80 1.01      626      730
##    y0[1]   52.88  52.86 1.14 1.11  50.99  54.74 1.01     1453     1062
##    y0[2]   30.94  30.93 1.13 1.10  29.08  32.77 1.00     1948     1884
##    y0[3]   57.29  57.29 1.19 1.16  55.33  59.23 1.01     1151     1367
##    y0[4]   40.25  40.21 1.09 1.02  38.39  42.05 1.00     1819     1556
##    y0[5]   57.12  57.13 1.21 1.16  55.11  58.97 1.00     1616     1731
##    y0[6]   36.53  36.53 1.10 1.10  34.68  38.29 1.00     1635     1466
##    y0[7]   26.82  26.81 1.16 1.13  24.91  28.73 1.02      359     1590
##    y0[8]   52.01  52.00 1.12 1.09  50.15  53.85 1.01     1183     1734
##    y0[9]   44.05  44.03 1.08 1.03  42.26  45.86 1.00     2022     1929
##    y0[10]  25.38  25.39 1.17 1.11  23.52  27.29 1.02     1109     1552
##    y0[11]  36.78  36.81 1.11 1.09  34.92  38.54 1.01     2022     1752
##    y0[12]  52.08  52.09 1.23 1.17  50.04  54.07 1.01      836      769
##    y0[13]  39.76  39.78 1.11 1.07  37.91  41.59 1.00     2077     1810
##    y0[14]  35.12  35.12 1.13 1.08  33.24  36.96 1.00     1878     1773
##    y0[15]  34.66  34.65 1.15 1.11  32.85  36.58 1.00     1958     1664
##    y0[16]  55.69  55.68 1.15 1.09  53.78  57.58 1.00     1407     1515
##    y0[17]  57.02  57.04 1.16 1.11  55.19  58.94 1.00     1593     1650
##    y0[18]  58.15  58.14 1.17 1.14  56.20  60.03 1.01     1417     1557
##    y0[19]  59.46  59.45 1.13 1.08  57.63  61.33 1.00     1761     1804
##    y0[20]  44.25  44.26 1.10 1.07  42.40  46.01 1.01     1296     1680
##    m1[1]   52.06  52.04 0.63 0.58  51.03  53.17 1.03      120      196
##    m1[2]   60.09  60.07 0.42 0.38  59.43  60.79 1.01     1170     1170
##    m1[3]   59.03  59.04 0.45 0.44  58.32  59.76 1.04       83     1009
##    m1[4]   54.49  54.48 0.45 0.45  53.76  55.23 1.05       62      518
##    m1[5]   48.89  48.89 0.38 0.37  48.29  49.52 1.03      104      511
##    m1[6]   43.26  43.26 0.31 0.29  42.77  43.78 1.00      758      828
##    m1[7]   38.02  38.03 0.30 0.29  37.51  38.52 1.01     1191      745
##    m1[8]   33.30  33.29 0.36 0.34  32.72  33.87 1.04       85      636
##    m1[9]   29.11  29.10 0.43 0.42  28.41  29.85 1.07       42      499
##    m1[10]  25.42  25.42 0.50 0.49  24.63  26.29 1.08       34      658
##    y1[1]   52.02  52.05 1.20 1.15  50.04  53.98 1.01     1214     1736
##    y1[2]   60.10  60.11 1.17 1.11  58.14  61.98 1.00     1774     1909
##    y1[3]   59.03  59.03 1.12 1.12  57.21  60.84 1.01     1417     1640
##    y1[4]   54.52  54.52 1.19 1.13  52.66  56.42 1.01     1313     1727
##    y1[5]   48.88  48.86 1.14 1.07  47.00  50.65 1.00     1703     1932
##    y1[6]   43.28  43.28 1.11 1.05  41.48  45.06 1.00     1886     1804
##    y1[7]   38.02  38.02 1.08 1.05  36.17  39.80 1.00     1948     1690
##    y1[8]   33.25  33.25 1.11 1.07  31.45  35.05 1.01     1678     1777
##    y1[9]   29.12  29.14 1.12 1.08  27.28  31.01 1.01     1496     1733
##    y1[10]  25.41  25.45 1.18 1.13  23.45  27.33 1.02      687     1133